Martinig, A. R., S. L. P. Burk, Y. Yang, M. Lagisz, and S. Nakagawa
Published
April 23, 2025
Things to do for April
Rerun all the models
Add some text so the reader can understand
Go through code and add more annotations
Ask questions about things you do not understand
1 Setting up
1.1 Load packages
Code
#force install pacman and orchaRd#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)#load all other packagespacman::p_load(tidyverse, # tidy family and related pacakges below kableExtra, # nice tables MuMIn, # multi-model inference emmeans, # post-hoc tests gridExtra, # may not use this pander, # nice tables metafor, # package for meta-analysis ape, # phylogenetic analysis ggtree, ggstance, MuMIn, # multi-model inference patchwork, # putting ggplots together - you need to install via devtool here, # making reading files easy orchaRd, # plotting orchard plots matrixcalc, # matrix calculations cowplot #multipanel plots)eval(metafor:::.MuMIn)
Table S1. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for intercept meta-analytic model with all five random effects.
Table S2. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for reduced intercept meta-analytic model without phylogeny as a random effect.
orchard_plot(mod_class, mod ="species_class", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S2. Mean fitness effect of different taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_class <-all_models(mod_class, mod ="species_class", type ="normal")# saving tab_class as RDSsaveRDS(tab_class, here("Rdata", "tab_class.RDS"))
Table S3. Mean estimates of different taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_type, mod ="study_type", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S3. Mean fitness effect of natural and semi-natural study types. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_type <-all_models(mod_type, mod ="study_type", type ="normal")# saving tab_type as RDSsaveRDS(tab_type, here("Rdata", "tab_type.RDS"))
Table S4. Mean estimates of natural and semi-natural study types with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_design, mod ="study_design", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S4. Mean fitness effect of gradient and contrast study designs. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_design <-all_models(mod_design, mod ="study_design", type ="normal")# saving tab_design as RDSsaveRDS(tab_design, here("Rdata", "tab_design.RDS"))
Table S5. Mean estimates of gradient and contrast study designs with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_disp, mod ="dispersal_type", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S5. Mean fitness effect of natal and breeding dispersal (or both dispersal types grouped). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_disp <-all_models(mod_disp, mod ="dispersal_type", type ="normal")# saving tab_disp as RDSsaveRDS(tab_disp, here("Rdata", "tab_disp.RDS"))
Table S6. Mean estimates of natal and breeding dispersal (or both dispersal types combined) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_phase, mod ="dispersal_phase", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S6. Mean fitness effect of prospection and settlement dispersal phases. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_phase <-all_models(mod_phase, mod ="dispersal_phase", type ="normal")# saving tab_phase as RDSsaveRDS(tab_phase, here("Rdata", "tab_phase.RDS"))
Table S7. Mean estimates of prospection and settlement dispersal phases with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_sex, mod ="sex", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S7. Mean fitness effect of females, males, or both sexes grouped. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex <-all_models(mod_sex, mod ="sex", type ="normal")# saving tab_sex as RDSsaveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))
Table S8. Mean estimates of females, males, or both sexes grouped with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_age1, mod ="age_class_clean", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S8. Mean fitness effect of life history stage (adult, juvenile, or mixed). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_age1 <-all_models(mod_age1, mod ="age_class_clean", type ="normal")# saving tab_age1 as RDSsaveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))
Table S9. Mean estimates of life history stage (adult, juvenile, or mixed) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_age2, mod ="age_class", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S9. Mean fitness effect of life history stage at the level of precision reported (A: adult, Y: yearling, J: juvenile). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)
3.3 Higher level
Code
mod_fit1 <-rma.mv(yi = yi, V = VCV,mod =~ fitness_higher_level -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_fit1)
orchard_plot(mod_fit1, mod ="fitness_higher_level", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S10. Mean fitness effect of different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_fit1 <-all_models(mod_fit1, mod ="fitness_higher_level", type ="normal")# saving tab_fit1 as RDSsaveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))
Table S10. Mean estimates of different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_fit2, mod ="fitness_metric_clean", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S11. Mean fitness effect of different fitness components at a finer resolution. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)
Code
mod_gen <-rma.mv(yi = yi, V = VCV,mod =~ whose_fitness -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_gen)
orchard_plot(mod_gen, mod ="whose_fitness", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)
Figure S12. Mean fitness effect of focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_gen <-all_models(mod_gen, mod ="whose_fitness", type ="normal")# saving tab_gen as RDSsaveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))
Table S11. Mean estimates of non-dispersers and dispersers and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
# getting absolute values for effect sizes as we expect shorter years will have more fluctuationsdat$ln_study_duration <-log(dat$study_duration+1)mod_dur <-rma.mv(yi = yi, V = VCV,mod =~ ln_study_duration,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_dur)
Figure S13. Mean fitness effect of study duration (log years). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Figure S14. Mean fitness effect of study duration (log years) from absolute value analysis. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
3.6 Normal analysis
Code
mod_focus <-rma.mv(yi = yi, V = VCV,mod =~ fitness_main_focus -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_focus)
orchard_plot(mod_focus, mod ="fitness_main_focus", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)
Figure S15. Mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_focus <-all_models(mod_focus, mod ="fitness_main_focus", type ="normal")# saving tab_focus as RDSsaveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))
Table S14. Mean estimates of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
# homoorchard_plot(mod_focus1, mod ="fitness_main_focus", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)
Figure S16. Maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# heteroorchard_plot(mod_focus2, mod ="fitness_main_focus", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)
Figure S17. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
3.8 Normal analysis
Code
mod_confirm <-rma.mv(yi = yi, V = VCV,mod =~ confirmation_bias -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_confirm)
orchard_plot(mod_confirm, mod ="confirmation_bias", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)
Figure S18. Mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_confirm <-all_models(mod_confirm, mod ="confirmation_bias", type ="normal")# saving tab_confirm as RDSsaveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))
Table S15. Mean estimates of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
# homoorchard_plot(mod_confirm1, mod ="confirmation_bias", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)
Figure S19. Maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# heteroorchard_plot(mod_confirm2, mod ="confirmation_bias", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)
Figure S20. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach dat$sex_species_class <-as.factor(paste(dat$sex, dat$species_class ,sep ="_"))mod_sex_species_class <-rma.mv(yi = yi, V = VCV,mod =~ sex_species_class -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_sex_species_class)
orchard_plot(mod_sex_species_class, mod ="sex_species_class", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S21. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_species_class <-all_models(mod_sex_species_class, mod ="sex_species_class", type ="normal")# saving tab_sex_species_class as RDSsaveRDS(tab_sex_species_class, here("Rdata", "tab_sex_species_class.RDS"))
Table S16. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_sex_whose_fitness, mod ="sex_whose_fitness", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S22. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_whose_fitness <-all_models(mod_sex_whose_fitness, mod ="sex_whose_fitness", type ="normal")# saving tab_sex_whose_fitness as RDSsaveRDS(tab_sex_whose_fitness, here("Rdata", "tab_sex_whose_fitness.RDS"))
Table S17. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_sex_fitness, mod ="sex_fitness_higher_level", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)
Figure S23. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_fitness <-all_models(mod_sex_fitness, mod ="sex_fitness_higher_level", type ="normal")# saving tab_sex_fitness as RDSsaveRDS(tab_sex_fitness, here("Rdata", "tab_sex_fitness.RDS"))
Table S18. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
orchard_plot(mod_sex_dispersal_type, mod ="sex_dispersal_type", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, trunk.size =0.5, angle =45)
Figure S24. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types combined). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_dispersal_type <-all_models(mod_sex_dispersal_type, mod ="sex_dispersal_type", type ="normal")# saving tab_sex_dispersal_type as RDSsaveRDS(tab_sex_dispersal_type, here("Rdata", "tab_sex_dispersal_type.RDS"))
Table S19. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types grouped) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).
There are two steps, which can be done using phylo_model().
5.1 Estimate species-level effect sizes.
Code
# function to estimate species-level effect sizesphylo_model <-function(data = dat, tree, correlation =0.5) {# accounting for within-species correlation Vm <-vcalc(vi = data$vi, cluster = data$shared_group, obs = data$effectID, rho = correlation) mod <-rma.mv(yi = yi, V = Vm,mod =~1,data = data,random =list(~1| effectID,~1| paperID,~1| phylogeny),#R= list(phylogeny = cor_tree),test ="t",sparse =TRUE)# get random effects# dat3 <- dat2[c("phylogeny", "yi", "vi")] REs <-ranef(mod)$phylogeny# calculate means and CIs dat <-data.frame(phylogeny =rownames(REs),estimate = mod$beta[1] + REs$intrcpt,se =sqrt((mod$se[1])^2+ (REs$se)^2)) %>%mutate(lb = estimate - se *qnorm(1-0.05/2), ub = estimate + se *qnorm(1-0.05/2)) %>%arrange(phylogeny)# get species class data species_class <- dplyr::distinct(data, phylogeny, .keep_all =TRUE) %>% dplyr::select(species_class, phylogeny)# extract tip labels from the tree and merge with species class tip.label <-data.frame(phylogeny = tree$tip.label) # extract tip label species_class2 <-left_join(tip.label, species_class, by ="phylogeny")return(list(data = dat, class = species_class2)) # return processed data and species class}# estimate species-level effect sizes# you can specify your own correlation coefficient via `correlation`. For illustration purpose, I used 0.5.result <-phylo_model(data = dat, tree = tree, correlation =0.5)# find the significant speciesresult$data <- result$data %>%mutate(p_value =2* (1-pnorm(abs(estimate / se))),significance =ifelse(p_value <0.05, "Significant", "NonSignificant") )result$data <- result$data %>%mutate(significance_label =ifelse(p_value <0.05, "*", NA))
5.2 Plot tree and add species-level effect size estimates
# adjust widths of the plot facets as necessary, if it improves visualization #| fig-cap: "**Figure S25.** Phylogenetic tree (left panel) and species-level effect size estimates (right panel). Significant relationships are indicated with a red asterisk. Dispersal was positive for ruffed grouse (Bonasa umbellus; Zr = 0.42, 95% CIs: 0.16 to 0.68) and negative for brown jays (Cyanocorax morio; Zr = -0.30, 95% CIs: -0.58 to -0.01), Tasmanian nativehens (Gallinula mortierii; Zr = -0.34, 95% CIs: -0.60 to X), and cougars (Puma concolor; Zr = -0.61, 95% CIs: -0.90 to -0.32)."facet_widths(p3, c(Tree =0.6, `Effect size`=0.4))
# relative importance of each predictor for all the modelsimportance <-sw(candidates)importance
age_class_clean dispersal_phase whose_fitness
Sum of weights: 0.74 0.40 0.40
N containing models: 32 32 32
fitness_higher_level sex dispersal_type
Sum of weights: 0.28 0.20 0.14
N containing models: 32 32 32
# if each group n is available - assume n/2dat$effectN <-ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n, (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))dat$sqeffectN <-sqrt(dat$effectN)mod_egger <-rma.mv(yi = yi, V = VCV,mod =~ sqeffectN,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_egger)
small <-bubble_plot(mod_egger,mod ="sqeffectN",group ="paperID",xlab ="sqrt(Effective N)",ylab ="Effect Size: Zr",g =TRUE)
Code
small
Figure S26. Relationship between sample size and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# creating publication year from "reference"dat$reference <-as.character(dat$reference)dat$year <-with(dat, substr(reference, nchar(reference)-4, nchar(reference)))dat$year <-as.integer(ifelse(dat$year =="2017a"| dat$year =="2017b", 2017, dat$year))# decline effectmod_dec <-rma.mv(yi = yi, V = vi,mod =~ year,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_dec)
Figure S27. Relationship between publication year and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# All togethermod_comb <-rma.mv(yi = yi, V = VCV,mod =~ year + sqeffectN,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),# ~ 1 | phylogeny),# R= list(phylogeny = cor_tree),test ="t",sparse =TRUE)summary(mod_comb)
Figure S28. Combined effect of sample size and publication year on effect size (plotting sample size only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Figure S29. Combined effect of sample size and publication year on effect size (plotting publication year only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
dat <- dat %>%mutate(leave_out = reference)dat$leave_out <-as.factor(dat$leave_out)LeaveOneOut_effectsize <-list()for (i in1:length(levels(dat$leave_out))) { temp_dat <- dat %>%filter(leave_out !=levels(dat$leave_out)[i]) VCV_leaveout <-vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho =0.5) LeaveOneOut_effectsize[[i]] <-rma.mv(yi = yi,V = VCV_leaveout, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE,data = temp_dat[temp_dat$leave_out !=levels(temp_dat$leave_out)[i], ])}# writing function for extracting est, ci.lb, and ci.ub from all modelsest.func <-function(model) { df <-data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)return(df)}# using dplyr to form data frameMA_oneout <-lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>% bind_rows %>%mutate(left_out =levels(dat$leave_out))# telling ggplot to stop reordering factorsMA_oneout$left_out <-as.factor(MA_oneout$left_out)MA_oneout$left_out <-factor(MA_oneout$left_out, levels = MA_oneout$left_out)# saving the runssaveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))
---title: "**Challenging assumptions: meta-analysis reveals dispersal costs are not universal.**"subtitle: "**Supplmentary Material**"author: "Martinig, A. R., S. L. P. Burk, Y. Yang, M. Lagisz, and S. Nakagawa"date: "`r Sys.Date()`"format: html: toc: true toc-location: left toc-depth: 3 toc-title: "**Table of Contents**" output-file: "index.html" theme: dark: darkly light: flatly embed-resources: true code-fold: true code-tools: true number-sections: true #bibliography: ./bib/ref.bib fontsize: "12" max-width: "10" code-overflow: wrapcrossref: fig-title: Figure # (default is "Figure") tbl-title: Table # (default is "Table") title-delim: — # (default is ":") fig-prefix: Fig. # (default is "Figure") tbl-prefix: Tab. # (default is "Table")editor_options: chunk_output_type: consoleexecute: warning: false message: false tidy: true #cache: true---Things to do for April1. Rerun all the models3. Add some text so the reader can understand4. Go through code and add more annotations5. Ask questions about things you do not understand## Setting up### Load packages```{r load packages}#force install pacman and orchaRd#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)#load all other packagespacman::p_load(tidyverse, # tidy family and related pacakges below kableExtra, # nice tables MuMIn, # multi-model inference emmeans, # post-hoc tests gridExtra, # may not use this pander, # nice tables metafor, # package for meta-analysis ape, # phylogenetic analysis ggtree, ggstance, MuMIn, # multi-model inference patchwork, # putting ggplots together - you need to install via devtool here, # making reading files easy orchaRd, # plotting orchard plots matrixcalc, # matrix calculations cowplot #multipanel plots)eval(metafor:::.MuMIn)```### Load data: paper and tree data```{r load data}#the dataset with effect sizesdat <- read.csv(here("data", "clean_data.csv"))#dim(dat)#head(dat)#the phylogenetic treetree <- read.tree(here("data", "species_tree.tre"))```### Calculate effect sizes and sampling variances using custom functions```{r calculate effect sizes}# function to calculate effect sizes# Zr - correlation# there is always neffect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments est , se, p_val, direction, method){ if(method == "mean_method"){ h <- n/n1 + n/n2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) ) j <- 1 - (3 / (4*n - 9)) d <- ((m2 - m1) / s_pool) * j r_pb <- d / sqrt(d^2 + h) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "count_method"){ h <- n/n1 + n/n2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 s1 <- sqrt(m1) s2 <- sqrt(m2) s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) ) j <- 1 - (3 / (4*n - 9)) d <- ((m2 - m1) / s_pool) * j r_pb <- d / sqrt(d^2 + h) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "percent_method1"){ h <- n/n1 + n/n2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 m1 <- m1/100 m2 <- m2/100 s1 <- 1/sqrt(8) s2 <- 1/sqrt(8) m1 <- asin(sqrt(m1/100)) m2 <- asin(sqrt(m2/100)) s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) ) j <- 1 - (3 / (4*n - 9)) d <- ((m2 - m1) / s_pool) * j r_pb <- d / sqrt(d^2 + h) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "percent_method2"){ h <- n/n1 + n/n2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 m1 <- m1/100 m2 <- m2/100 sd1 <- sd1/100 sd2 <- sd2/100 s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1))) s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2))) m1 <- asin(sqrt(m1/100)) m2 <- asin(sqrt(m2/100)) s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) ) j <- 1 - (3 / (4*n - 9)) d <- ((m2 - m1) / s_pool) * j r_pb <- d / sqrt(d^2 + h) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "proportion_method1"){ h <- n/n1 + n/n2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 s1 <- 1/sqrt(8) s2 <- 1/sqrt(8) m1 <- asin(sqrt(m1)) m2 <- asin(sqrt(m2)) s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) ) j <- 1 - (3 / (4*n - 9)) d <- ((m2 - m1) / s_pool) * j r_pb <- d / sqrt(d^2 + h) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "proportion_method2"){ h <- n/n1 + n/n2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1))) s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2))) m1 <- asin(sqrt(m1/100)) m2 <- asin(sqrt(m2/100)) s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) ) j <- 1 - (3 / (4*n - 9)) d <- ((m2 - m1) / s_pool) * j r_pb <- d / sqrt(d^2 + h) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "t_method1"){ p <- n1/n # prop for n1 q <- n2/n # prop for n2 r_pb <- est/sqrt(est^2 + n - 2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "t_method2"){ n1 <- n/2 n2 <- n/2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 r_pb <- est/sqrt(est^2 + n - 2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "F_method1"){ p <- n1/n # prop for n1 q <- n2/n # prop for n2 r_pb <- sqrt(est)/sqrt(est + n -2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r_b = r_b*(direction) r <- r_b }else if(method == "F_method2"){ n1 <- n/2 n2 <- n/2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 r_pb <- sqrt(est)/sqrt(est + n -2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r_b = r_b*(direction) r <- r_b }else if(method == "p_method1"){ p <- n1/n # prop for n1 q <- n2/n # prop for n2 t <- qt(1 - p_val, n - 2) r_pb <- t/sqrt(t^2 + n -2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r_b <- r_b*(direction) r <- r_b }else if(method == "p_method2"){ n1 <- n/2 n2 <- n/2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 t <- qt(1 - p_val, n - 2) r_pb <- t/sqrt(t^2 + n -2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r_b <- r_b*(direction) r <- r_b }else if(method == "correlation_method1"){ r <- est }else if(method == "correlation_method2"){ r <- 2*sin((pi/6)*est) }else if(method == "estimate_method1"){ p <- n1/n # prop for n1 q <- n2/n # prop for n2 t <- est/se r_pb <- t/sqrt(t^2+ n -2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r }else if(method == "estimate_method2"){ n1 <- n/2 n2 <- n/2 p <- n1/n # prop for n1 q <- n2/n # prop for n2 t <- est/se r_pb <- t/sqrt(t^2+ n -2) r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p))) r <- r_b #r_b = r } if(r >= 1){ # if over 1, we use 0.99 Zr <- atanh(0.99) }else if(r <= -1){ Zr <- atanh(-0.99) # r = correlation } else { Zr <- atanh(r) # r = correlation } VZr <- 1 /(n - 3) data.frame(ri = r, yi = Zr , vi = VZr)}###########' Title: Contrast name generator#'#' @param name: a vector of character stringscont_gen <- function(name) { combination <- combn(name, 2) name_dat <- t(combination) names <- paste(name_dat[, 1], name_dat[, 2], sep = "-") return(names)}#' @title get_pred1: intercept-less model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred1 <- function (model, mod = " ") { name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, ""))) len <- length(name) if (len != 1) { newdata <- diag(len) pred <- metafor::predict.rma(model, newmods = newdata, tau2.levels = 1:len) } else { pred <- metafor::predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate, lowerCL = lowerCL, upperCL = upperCL, pval = model$pval, lowerPR = lowerPR, upperPR = upperPR)}#' @title get_pred2: normal model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred2 <- function (model, mod = " ") { name <- as.factor(str_replace(row.names(model$beta), paste0("relevel", "\\(", mod,", ref = name","\\)"),"")) len <- length(name) if(len != 1){ newdata <- diag(len) pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1]) } else { pred <- predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate, lowerCL = lowerCL, upperCL = upperCL, pval = model$pval, lowerPR = lowerPR, upperPR = upperPR)}#' @title mr_results#' @description Function to put results of meta-regression and its contrasts#' @param res1: data frame 1#' @param res1: data frame 2mr_results <- function(res1, res2) { restuls <-tibble( `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)), Estimate = c(res1$estimate, res2$estimate), `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL), `Upper CI [0.975]` = c(res1$upperCL, res2$upperCL), `P value` = c(res1$pval, res2$pval), `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR), `Upper PI [0.975]` = c(res1$upperPR, res2$upperPR), )}#' @title all_models#' @description Function to take all possible models and get their results#' @param model: intercept-less model#' @param mod: the name of a moderator all_models <- function(model, mod = " ", type = "normal") { # getting the level names out level_names <- levels(factor(model$data[[mod]])) dat2 <- model$data mod <- mod # creating a function to run models run_rma1 <- function(name) { # variance covarinace matrix for sampling errors VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5) #VCV1 <- nearPD(VCV1)$mat rma.mv(yi, V = VCV1, mods = ~relevel(dat2[[mod]], ref = name), random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), data = dat2, test = "t", sparse = TRUE, control = list(optimizer = "Nelder-Mead")) } run_rma2 <- function(name) { # variance covarinace matrix for sampling errors VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5) #VCVa <- nearPD(VCVa)$mat rma.mv(abs_yi2, V = VCVa, mods = ~relevel(dat2[[mod]], ref = name), random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), data = dat2, test = "t", sparse = TRUE, control = list(optimizer = "Nelder-Mead")) } run_rma3 <- function(name) { # variance covarinace matrix for sampling errors VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5) #VCV1 <- nearPD(VCV1)$mat rma.mv(yi, V = VCV1, mods = ~relevel(dat2[[mod]], ref = name), random = list( # hetero in relation to mod formula(paste("~", mod, "| effectID")), ~ 1 | paperID, ~ 1 | species_cleaned), struct = "DIAG", data = dat2, test = "t", sparse = TRUE, control = list(optimizer = "Nelder-Mead")) }# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])# use different specifications for aboslute and hetero models if (type == "normal"){ model_all <- purrr::map(level_names[-length(level_names)], run_rma1) }else if(type == "absolute"){ model_all <- purrr::map(level_names[-length(level_names)], run_rma2) }else if(type == "hetero"){ model_all <- purrr::map(level_names[-length(level_names)], run_rma3) } # getting estimates from intercept-less models (means for all the groups) res1 <- get_pred1(model, mod = mod) # getting estiamtes from all contrast models res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod)) # a list of the numbers to take out unnecessary contrasts contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1)) res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) # creating a table res_tab <- mr_results(res1, res2) %>% kable("html", digits = 3) %>% kable_styling("striped", position = "left") %>% scroll_box(width = "100%") # results res_tab}# absolute value analyses# folded meanfolded_mu <-function(mean, variance){ mu <- mean sigma <- sqrt(variance) fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma)) fold_mu} # folded variancefolded_v <-function(mean, variance){ mu <- mean sigma <- sqrt(variance) fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma)) fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2) # adding se to make bigger mean fold_v <-fold_se^2 fold_v} ```### Preparing final dataset```{r final dataset}#| include: false#| # adding effect sizes to dataseteffect2 <- pmap_dfr(list(dat$mean_group_1, dat$mean_group_2, dat$variance_group_1, dat$variance_group_2, dat$n_group_1, dat$n_group_2, dat$n, dat$effect_size_value, dat$effect_size_variance, dat$effect_size_p_value_numeric, dat$direction_change, factor(dat$function_needed)), effect_size) # merging two data framesdat <- cbind(dat, effect2)# renaming X to effectIDcolnames(dat)[colnames(dat) == "X"] <- "effectID"# creating the phylogeny columndat$phylogeny <- gsub(" ", "_", dat$species_cleaned)# absolute values for the effect sizesdat$yi_abs <- folded_mu(dat$yi, dat$vi)dat$vi_abs <- folded_v(dat$yi, dat$vi)# checking species names match between tree and dataset# match(unique(dat$phylogeny), tree$tip.label)# match(tree$tip.label, unique(dat$phylogeny))# # intersect(unique(dat$phylogeny), tree$tip.label)# setdiff(unique(dat$phylogeny), tree$tip.label)# # # match(unique(dat$phylogeny), tree$tip.label)# sum(is.na(match(unique(dat$phylogeny), tree$tip.label)))# creating a correlation matrix from a treetree <- compute.brlen(tree)cor_tree <- vcv(tree,corr=T) # creating a VCV for sampling error variances # we assure shared groups will have r = 0.5# it is fine if these are not positive definite but changing thisVCV <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)#is.positive.definite(as.matrix(VCV), tol=1e-8)# focing the matrix to be positive definite# this should be in the method section#VCV <- nearPD(VCV)$mat# VCV for aboslute valuesVCVa <- vcalc(vi = dat$vi_abs, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)# marking all character strings in dat into factordat <- dat %>% mutate(across(where(is.character), as.factor))```## Intercept meta-analytic model### Phylogenetic multilevel models::: panel-tabset#### With all five random effects```{r phylogenetic multilevel models}#| eval: false#| # takes a while to runmod <- rma.mv(yi = yi, V = VCV, mod = ~ 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned, ~ 1 | phylogeny), R= list(phylogeny = cor_tree), test = "t", sparse = TRUE)# saving the runssaveRDS(mod, here("Rdata", "mod.RDS"))``````{r}# getting the saved modelmod <-readRDS(here("Rdata", "mod.RDS"))summary(mod)# Phylogeney accounts for nothing (0.00%), so we can take it out. round(i2_ml(mod),2) ``````{r}pred_mod <-predict.rma(mod) estimate <- pred_mod$pred lowerCL <- pred_mod$ci.lb upperCL <- pred_mod$ci.ub lowerPR <- pred_mod$cr.lb upperPR <- pred_mod$cr.ub table_mod <-tibble("Fixed effect"="intercept", estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = mod$pval,lowerPR = lowerPR, upperPR = upperPR)# creating a tabletab_mod <- table_mod %>%kable("html", digits =3) %>%kable_styling("striped", position ="left") %>%scroll_box(width ="100%")# saving tab_mod as RDSsaveRDS(tab_mod, here("Rdata", "tab_mod.RDS"))```**Table S1.** Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for intercept meta-analytic model with all five random effects.```{r}tab_mod <-readRDS(here("Rdata", "tab_mod.RDS"))tab_mod```#### Reduced multilevel modelWe remove phylogeny as a random effect from our meta-analytic model because phylogeny accounts for nothing (0%)```{r}mod1 <-rma.mv(yi = yi, V = VCV,mod =~1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod1)round(i2_ml(mod1),2) ``````{r}#| fig-cap: "**Figure S1.** Overall meta-analytic mean effect size (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE)."reduced_intercept<-orchard_plot(mod1, xlab ="Effect Size: Zr", group ="paperID", branch.size =4)``````{r}pred_mod1 <-predict.rma(mod1) estimate <- pred_mod1$pred lowerCL <- pred_mod1$ci.lb upperCL <- pred_mod1$ci.ub lowerPR <- pred_mod1$cr.lb upperPR <- pred_mod1$cr.ub table_mod1 <-tibble("Fixed effect"="intercept", estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = mod1$pval,lowerPR = lowerPR, upperPR = upperPR)# creating a tabletab_mod1 <- table_mod1 %>%kable("html", digits =3) %>%kable_styling("striped", position ="left") %>%scroll_box(width ="100%")# saving tab_mod1 as RDSsaveRDS(tab_mod1, here("Rdata", "tab_mod1.RDS"))```**Table S2.** Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for reduced intercept meta-analytic model without phylogeny as a random effect.```{r}tab_mod1 <-readRDS(here("Rdata", "tab_mod1.RDS"))tab_mod1```:::## Uni-moderator meta-regression models For each of our moderators, we run a uni-moderator meta-regression model::: panel-tabset## Species class```{r species class}mod_class <- rma.mv(yi = yi, V = VCV, mod = ~ species_class - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_class)r2_ml(mod_class)``````{r}#| fig-cap: "**Figure S2.** Mean fitness effect of different taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_class, mod ="species_class", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_class <-all_models(mod_class, mod ="species_class", type ="normal")# saving tab_class as RDSsaveRDS(tab_class, here("Rdata", "tab_class.RDS"))```**Table S3.** Mean estimates of different taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_class <-readRDS(here("Rdata", "tab_class.RDS"))tab_class```## Study type```{r study type}mod_type <- rma.mv(yi = yi, V = VCV, mod = ~ study_type - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_type)r2_ml(mod_type)``````{r}#| fig-cap: "**Figure S3.** Mean fitness effect of natural and semi-natural study types. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_type, mod ="study_type", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)``````{r}#| eval: falsetab_type <-all_models(mod_type, mod ="study_type", type ="normal")# saving tab_type as RDSsaveRDS(tab_type, here("Rdata", "tab_type.RDS"))```**Table S4.** Mean estimates of natural and semi-natural study types with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_type <-readRDS(here("Rdata", "tab_type.RDS"))tab_type```## Study design```{r study design}mod_design <- rma.mv(yi = yi, V = VCV, mod = ~ study_design - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_design)r2_ml(mod_design)``````{r}#| fig-cap: "**Figure S4.** Mean fitness effect of gradient and contrast study designs. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_design, mod ="study_design", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_design <-all_models(mod_design, mod ="study_design", type ="normal")# saving tab_design as RDSsaveRDS(tab_design, here("Rdata", "tab_design.RDS"))```**Table S5.** Mean estimates of gradient and contrast study designs with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_design <-readRDS(here("Rdata", "tab_design.RDS"))tab_design```## Dispersal type```{r dispersal type}mod_disp <- rma.mv(yi = yi, V = VCV, mod = ~ dispersal_type - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_disp)r2_ml(mod_disp)``````{r}#| fig-cap: "**Figure S5.** Mean fitness effect of natal and breeding dispersal (or both dispersal types grouped). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_disp, mod ="dispersal_type", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_disp <-all_models(mod_disp, mod ="dispersal_type", type ="normal")# saving tab_disp as RDSsaveRDS(tab_disp, here("Rdata", "tab_disp.RDS"))```**Table S6.** Mean estimates of natal and breeding dispersal (or both dispersal types combined) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_disp <-readRDS(here("Rdata", "tab_disp.RDS"))tab_disp```## Dispersal phase```{r dispersal phase}mod_phase <- rma.mv(yi = yi, V = VCV, mod = ~ dispersal_phase - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_phase)r2_ml(mod_phase)``````{r}#| fig-cap: "**Figure S6.** Mean fitness effect of prospection and settlement dispersal phases. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_phase, mod ="dispersal_phase", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_phase <-all_models(mod_phase, mod ="dispersal_phase", type ="normal")# saving tab_phase as RDSsaveRDS(tab_phase, here("Rdata", "tab_phase.RDS"))```**Table S7.** Mean estimates of prospection and settlement dispersal phases with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_phase <-readRDS(here("Rdata", "tab_phase.RDS"))tab_phase```## Sex```{r sex}mod_sex <- rma.mv(yi = yi, V = VCV, mod = ~ sex - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_sex)r2_ml(mod_sex)``````{r}#| fig-cap: "**Figure S7.** Mean fitness effect of females, males, or both sexes grouped. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_sex, mod ="sex", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)``````{r}#| eval: falsetab_sex <-all_models(mod_sex, mod ="sex", type ="normal")# saving tab_sex as RDSsaveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))```**Table S8.** Mean estimates of females, males, or both sexes grouped with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))tab_sex```## Life history stage### Higher life history grouping```{r life history stage}mod_age1 <- rma.mv(yi = yi, V = VCV, mod = ~ age_class_clean - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_age1)r2_ml(mod_age1)``````{r}#| fig-cap: "**Figure S8.** Mean fitness effect of life history stage (adult, juvenile, or mixed). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_age1, mod ="age_class_clean", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)``````{r}#| eval: falsetab_age1 <-all_models(mod_age1, mod ="age_class_clean", type ="normal")# saving tab_age1 as RDSsaveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))```**Table S9.** Mean estimates of life history stage (adult, juvenile, or mixed) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))tab_age1```### Finer life stage grouping```{r}# age_classmod_age2 <-rma.mv(yi = yi,V = VCV,mod =~ age_class -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_age2)r2_ml(mod_age2)``````{r}#| fig-cap: "**Figure S9.** Mean fitness effect of life history stage at the level of precision reported (A: adult, Y: yearling, J: juvenile). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_age2, mod ="age_class", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)# table not created for this (not meaningful for some levels)```## Fitness metric### Higher level```{r fitness metric}mod_fit1 <- rma.mv(yi = yi, V = VCV, mod = ~ fitness_higher_level - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_fit1)r2_ml(mod_fit1)``````{r}#| fig-cap: "**Figure S10.** Mean fitness effect of different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_fit1, mod ="fitness_higher_level", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)``````{r}#| eval: falsetab_fit1 <-all_models(mod_fit1, mod ="fitness_higher_level", type ="normal")# saving tab_fit1 as RDSsaveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))```**Table S10.** Mean estimates of different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))tab_fit1```### Finer level```{r}mod_fit2 <-rma.mv(yi = yi, V = VCV,mod =~ fitness_metric_clean -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_fit2)r2_ml(mod_fit2)``````{r}#| fig-cap: "**Figure S11.** Mean fitness effect of different fitness components at a finer resolution. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_fit2, mod ="fitness_metric_clean", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)# table not created for this (not meaningful for some levels)```## Multi-generation```{r multi-generation}mod_gen <- rma.mv(yi = yi, V = VCV, mod = ~ whose_fitness - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_gen)r2_ml(mod_gen)``````{r}#| fig-cap: "**Figure S12.** Mean fitness effect of focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_gen, mod ="whose_fitness", xlab ="Effect Size: Zr", group ="paperID", branch.size =4)``````{r}#| eval: falsetab_gen <-all_models(mod_gen, mod ="whose_fitness", type ="normal")# saving tab_gen as RDSsaveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))```**Table S11.** Mean estimates of non-dispersers and dispersers and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))tab_gen```## Study duration### Normal analysis```{r study duration}# getting absolute values for effect sizes as we expect shorter years will have more fluctuationsdat$ln_study_duration <- log(dat$study_duration+1)mod_dur <- rma.mv(yi = yi, V = VCV, mod = ~ ln_study_duration, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_dur)r2_ml(mod_dur)``````{r}#| fig-cap: "**Figure S13.** Mean fitness effect of study duration (log years). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."bubble_plot(mod_dur,mod ="ln_study_duration",group ="paperID",xlab ="Study duration",ylab ="Effect Size: Zr",g =TRUE)``````{r}# absolute value analysismod_dur_abs <-rma.mv(yi = yi_abs,V = VCVa,mod =~ ln_study_duration,data = dat,random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",sparse =TRUE)summary(mod_dur_abs)r2_ml(mod_dur_abs)``````{r}#| fig-cap: "**Figure S14.** Mean fitness effect of study duration (log years) from absolute value analysis. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."bubble_plot(mod_dur_abs,mod ="ln_study_duration",group ="paperID",xlab ="Study duration",ylab ="Effect Size: Zr",g =TRUE)```## Selection bias### Normal analysis```{r selection bias}mod_focus <- rma.mv(yi = yi, V = VCV, mod = ~ fitness_main_focus - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_focus)r2_ml(mod_focus)``````{r}#| fig-cap: "**Figure S15.** Mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_focus, mod ="fitness_main_focus", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)``````{r}#| eval: falsetab_focus <-all_models(mod_focus, mod ="fitness_main_focus", type ="normal")# saving tab_focus as RDSsaveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))```**Table S14.** Mean estimates of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))tab_focus```### Comparing normal model to heteroscedastic model```{r}# homomod_focus1 <-rma.mv(yi = yi, V = VCV,mod =~ fitness_main_focus -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",method ="ML",sparse =TRUE)summary(mod_focus1)r2_ml(mod_focus1)# heteromod_focus2 <-rma.mv(yi = yi, V = VCV,mod =~ fitness_main_focus,data = dat, random =list(~ fitness_main_focus | effectID,~1| paperID,~1| species_cleaned),struct ="DIAG",method ="ML",test ="t",sparse =TRUE,control=list(optimizer="nlminb", rel.tol=1e-8))summary(mod_focus2)#r2_ml(mod_focusA1)# they are not significantly differentanova(mod_focus1, mod_focus2)``````{r}#| fig-cap: "**Figure S16.** Maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."# homoorchard_plot(mod_focus1, mod ="fitness_main_focus", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)``````{r}#| fig-cap: "**Figure S17.** Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."# heteroorchard_plot(mod_focus2, mod ="fitness_main_focus", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)```## Confirmation bias### Normal analysis```{r confirmation bias}mod_confirm <- rma.mv(yi = yi, V = VCV, mod = ~ confirmation_bias - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_confirm)r2_ml(mod_confirm)``````{r}#| fig-cap: "**Figure S18.** Mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_confirm, mod ="confirmation_bias", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)``````{r}#| eval: falsetab_confirm <-all_models(mod_confirm, mod ="confirmation_bias", type ="normal")# saving tab_confirm as RDSsaveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))```**Table S15.** Mean estimates of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_confirm <-readRDS(here("Rdata", "tab_confirm.RDS"))tab_confirm```### Comparing normal model to heteroscedastic model```{r}# homomod_confirm1 <-rma.mv(yi = yi, V = VCV,mod =~ confirmation_bias -1,data = dat, random =list(~1| effectID,~1| paperID,~1| species_cleaned),test ="t",method ="ML",sparse =TRUE)summary(mod_confirm1)r2_ml(mod_confirm1)# heteromod_confirm2 <-rma.mv(yi = yi, V = VCV,mod =~ confirmation_bias,data = dat, random =list(~ confirmation_bias | effectID,~1| paperID,~1| species_cleaned),struct ="DIAG",method ="ML",test ="t",sparse =TRUE,control=list(optimizer="nlminb", rel.tol=1e-8))summary(mod_confirm2)# they are not significantly differentanova(mod_confirm1, mod_confirm2)``````{r}#| fig-cap: "**Figure S19.** Maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."# homoorchard_plot(mod_confirm1, mod ="confirmation_bias", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)``````{r}#| fig-cap: "**Figure S20.** Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."# heteroorchard_plot(mod_confirm2, mod ="confirmation_bias", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =90)```:::## Interaction meta-regression models (2 moderators)::: panel-tabset## Sex x Species class```{r sex x species class}# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach dat$sex_species_class <- as.factor(paste(dat$sex, dat$species_class ,sep = "_"))mod_sex_species_class <- rma.mv(yi = yi, V = VCV, mod = ~ sex_species_class - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_sex_species_class)r2_ml(mod_sex_species_class)``````{r}#| fig-cap: "**Figure S21.** Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_sex_species_class, mod ="sex_species_class", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_sex_species_class <-all_models(mod_sex_species_class, mod ="sex_species_class", type ="normal")# saving tab_sex_species_class as RDSsaveRDS(tab_sex_species_class, here("Rdata", "tab_sex_species_class.RDS"))```**Table S16.** Mean estimates for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_sex_species_class <-readRDS(here("Rdata", "tab_sex_species_class.RDS"))tab_sex_species_class```## Sex x Whose fitness```{r sex x whose fitness}dat$sex_whose_fitness <- as.factor(paste(dat$sex, dat$whose_fitness ,sep = "_"))mod_sex_whose_fitness <- rma.mv(yi = yi, V = VCV, mod = ~ sex_whose_fitness - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_sex_whose_fitness)r2_ml(mod_sex_whose_fitness)``````{r}#| fig-cap: "**Figure S22.** Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_sex_whose_fitness, mod ="sex_whose_fitness", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_sex_whose_fitness <-all_models(mod_sex_whose_fitness, mod ="sex_whose_fitness", type ="normal")# saving tab_sex_whose_fitness as RDSsaveRDS(tab_sex_whose_fitness, here("Rdata", "tab_sex_whose_fitness.RDS"))```**Table S17.** Mean estimates for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_sex_whose_fitness <-readRDS(here("Rdata", "tab_sex_whose_fitness.RDS"))tab_sex_whose_fitness```## Sex x Fitness (higher level)```{r sex x fitness}dat$sex_fitness_higher_level <- as.factor(paste(dat$sex, dat$fitness_higher_level ,sep = "_"))mod_sex_fitness <- rma.mv(yi = yi, V = VCV, mod = ~ sex_fitness_higher_level - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_sex_fitness)r2_ml(mod_sex_fitness)``````{r}#| fig-cap: "**Figure S23.** Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_sex_fitness, mod ="sex_fitness_higher_level", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, angle =45)``````{r}#| eval: falsetab_sex_fitness <-all_models(mod_sex_fitness, mod ="sex_fitness_higher_level", type ="normal")# saving tab_sex_fitness as RDSsaveRDS(tab_sex_fitness, here("Rdata", "tab_sex_fitness.RDS"))```**Table S18.** Mean estimates for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_sex_fitness <-readRDS(here("Rdata", "tab_sex_fitness.RDS"))tab_sex_fitness```## Sex x Dispersal type```{r sex x dispersal type}dat$sex_dispersal_type <- as.factor(paste(dat$sex, dat$dispersal_type ,sep = "_"))mod_sex_dispersal_type <- rma.mv(yi = yi, V = VCV, mod = ~ sex_dispersal_type - 1, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_sex_dispersal_type)r2_ml(mod_sex_dispersal_type)``````{r}#| fig-cap: "**Figure S24.** Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types combined). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."orchard_plot(mod_sex_dispersal_type, mod ="sex_dispersal_type", xlab ="Effect Size: Zr", group ="paperID", branch.size =4, trunk.size =0.5, angle =45)``````{r}#| eval: falsetab_sex_dispersal_type <-all_models(mod_sex_dispersal_type, mod ="sex_dispersal_type", type ="normal")# saving tab_sex_dispersal_type as RDSsaveRDS(tab_sex_dispersal_type, here("Rdata", "tab_sex_dispersal_type.RDS"))```**Table S19.** Mean estimates for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types grouped) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).```{r}tab_sex_dispersal_type <-readRDS(here("Rdata", "tab_sex_dispersal_type.RDS"))tab_sex_dispersal_type```:::## Species-level modelsThere are two steps, which can be done using phylo_model().### Estimate species-level effect sizes.```{r species-level effect sizes}# function to estimate species-level effect sizesphylo_model <- function(data = dat, tree, correlation = 0.5) { # accounting for within-species correlation Vm <- vcalc(vi = data$vi, cluster = data$shared_group, obs = data$effectID, rho = correlation) mod <- rma.mv(yi = yi, V = Vm, mod = ~ 1, data = data, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | phylogeny), #R= list(phylogeny = cor_tree), test = "t", sparse = TRUE) # get random effects # dat3 <- dat2[c("phylogeny", "yi", "vi")] REs <- ranef(mod)$phylogeny # calculate means and CIs dat <- data.frame(phylogeny = rownames(REs), estimate = mod$beta[1] + REs$intrcpt, se = sqrt((mod$se[1])^2 + (REs$se)^2)) %>% mutate(lb = estimate - se * qnorm(1 - 0.05/2), ub = estimate + se * qnorm(1 - 0.05/2)) %>% arrange(phylogeny) # get species class data species_class <- dplyr::distinct(data, phylogeny, .keep_all = TRUE) %>% dplyr::select(species_class, phylogeny) # extract tip labels from the tree and merge with species class tip.label <- data.frame(phylogeny = tree$tip.label) # extract tip label species_class2 <- left_join(tip.label, species_class, by = "phylogeny") return(list(data = dat, class = species_class2)) # return processed data and species class}# estimate species-level effect sizes# you can specify your own correlation coefficient via `correlation`. For illustration purpose, I used 0.5.result <- phylo_model(data = dat, tree = tree, correlation = 0.5)# find the significant speciesresult$data <- result$data %>% mutate( p_value = 2 * (1 - pnorm(abs(estimate / se))), significance = ifelse(p_value < 0.05, "Significant", "NonSignificant") )result$data <- result$data %>% mutate(significance_label = ifelse(p_value < 0.05, "*", NA))```### Plot tree and add species-level effect size estimates```{r plot tree with species}# plot treep1 <- ggtree(tree, layout = "rectangular", cex = 0.3)p2 <- p1 %<+% result$class + geom_tiplab(aes(color = species_class), size = 3, align = T, offset = 0.05) + geom_tippoint(aes(color = species_class)) + guides(color = "none")p3 <- p2 + geom_facet( panel = "Effect size", data = result$data, geom = ggstance::geom_pointrangeh, mapping = aes(x = estimate, xmin = lb, xmax = ub, color = species_class) ) + geom_facet( panel = "Effect size", data = result$data, geom = geom_text, mapping = aes(x = lb, label = significance_label), hjust = 1.5, size = 3, color = "red" # adjust position of star * ) + geom_facet( panel = "Effect size", data = result$data, geom = geom_vline, mapping = aes(xintercept = 0), linetype = "dashed", color = "gray" ) + theme_tree2()``````{r, fig.width=20, fig.height=30}# adjust widths of the plot facets as necessary, if it improves visualization #| fig-cap: "**Figure S25.** Phylogenetic tree (left panel) and species-level effect size estimates (right panel). Significant relationships are indicated with a red asterisk. Dispersal was positive for ruffed grouse (Bonasa umbellus; Zr = 0.42, 95% CIs: 0.16 to 0.68) and negative for brown jays (Cyanocorax morio; Zr = -0.30, 95% CIs: -0.58 to -0.01), Tasmanian nativehens (Gallinula mortierii; Zr = -0.34, 95% CIs: -0.60 to X), and cougars (Puma concolor; Zr = -0.61, 95% CIs: -0.90 to -0.32)."facet_widths(p3, c(Tree = 0.6, `Effect size` = 0.4))```:::## Multi-moderator-regression model::: panel-tabset## Full model```{r multi-moderator-regression model}# used "maximum likelihood" for model selection######################## Multi-variable models#######################mod_full <- rma.mv(yi = yi, V = VCV, mod = ~ 1 + sex + age_class_clean + whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), method = "ML", test = "t", sparse = TRUE)summary(mod_full)round(r2_ml(mod_full)*100, 2)```## Multi-model inference```{r}#| eval: false# multi-model selectioncandidates <-dredge(mod_full, trace =2)# saving candidates as RDSsaveRDS(candidates, here("Rdata", "candidates.RDS"))``````{r}candidates <-readRDS(here("Rdata", "candidates.RDS"))candidates# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <2) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <2)) mr_averaged_aic2# relative importance of each predictor for all the modelsimportance <-sw(candidates)importance```:::## Publication bias::: panel-tabset## Egger regression: uni-moderator```{r Egger regression}# if each group n is available - assume n/2dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n, (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))dat$sqeffectN <- sqrt(dat$effectN)mod_egger <- rma.mv(yi = yi, V = VCV, mod = ~ sqeffectN, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_egger)round(r2_ml(mod_egger)*100, 2)``````{r}small <-bubble_plot(mod_egger,mod ="sqeffectN",group ="paperID",xlab ="sqrt(Effective N)",ylab ="Effect Size: Zr",g =TRUE)``````{r}#| fig-cap: "**Figure S26.** Relationship between sample size and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."small```## Decline effect: uni-moderator```{r decline effect}# creating publication year from "reference"dat$reference <- as.character(dat$reference)dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))# decline effectmod_dec <- rma.mv(yi = yi, V = vi, mod = ~ year, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE)summary(mod_dec)round(r2_ml(mod_dec)*100, 2)decline <- bubble_plot(mod_dec, mod = "year", group = "paperID", xlab = "Publication year", ylab = "Effect Size: Zr", g = TRUE)``````{r}#| fig-cap: "**Figure S27.** Relationship between publication year and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."decline ```## Combined analysis```{r combined analysis}# All togethermod_comb <- rma.mv(yi = yi, V = VCV, mod = ~ year + sqeffectN, data = dat, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), # ~ 1 | phylogeny), # R= list(phylogeny = cor_tree), test = "t", sparse = TRUE)summary(mod_comb)round(r2_ml(mod_comb)*100, 2)``````{r}#| fig-cap: "**Figure S28.** Combined effect of sample size and publication year on effect size (plotting sample size only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."# combined sample size and yearbubble_plot(mod_comb,mod ="sqeffectN",group ="paperID",xlab ="sqrt(Effective N)",ylab ="Effect Size: Zr",g =TRUE)``````{r}#| fig-cap: "**Figure S29.** Combined effect of sample size and publication year on effect size (plotting publication year only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies)."# combined sample size and yearbubble_plot(mod_comb,mod ="year",group ="paperID",xlab ="Publication year",ylab ="Effect Size: Zr",g =TRUE)```## Leave one study out (sensitivity analysis)```{r leave-1study-out}#| eval: falsedat <- dat %>% mutate(leave_out = reference)dat$leave_out <- as.factor(dat$leave_out)LeaveOneOut_effectsize <- list()for (i in 1:length(levels(dat$leave_out))) { temp_dat <- dat %>% filter(leave_out != levels(dat$leave_out)[i]) VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5) LeaveOneOut_effectsize[[i]] <- rma.mv(yi = yi, V = VCV_leaveout, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE, data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])}# writing function for extracting est, ci.lb, and ci.ub from all modelsest.func <- function(model) { df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub) return(df)}# using dplyr to form data frameMA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>% bind_rows %>% mutate(left_out = levels(dat$leave_out))# telling ggplot to stop reordering factorsMA_oneout$left_out <- as.factor(MA_oneout$left_out)MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)# saving the runssaveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))``````{r}#| fig-height: 16MA_oneout <-readRDS(here("Rdata", "MA_oneout.RDS"))# plottingleaveoneout <-ggplot(MA_oneout) +geom_hline(yintercept =0, lty =2, lwd =1) +geom_hline(yintercept = mod1$ci.lb, lty =3, lwd =0.75, colour ="black") +geom_hline(yintercept = mod1$b, lty =1, lwd =0.75, colour ="black") +geom_hline(yintercept = mod1$ci.ub,lty =3, lwd =0.75, colour ="black") +geom_pointrange(aes(x = left_out, y = est,ymin = lower, ymax = upper)) +xlab("Study left out") +ylab("Zr (effect size), 95% CI") +coord_flip() +theme(panel.grid.minor =element_blank()) +theme_bw() +theme(panel.grid.major =element_blank()) +theme(panel.grid.minor.x =element_blank()) +theme(axis.text.y =element_text(size =6))``````{r, fig.width=20, fig.height=30}#| fig-cap: "**Figure S30.** Leave one study out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals."leaveoneout```## Leave one species out (sensitivity analysis)```{r leave-1species-out}#| eval: falsedat <- dat %>% mutate(leave_out = species_cleaned)dat$leave_out <- as.factor(dat$leave_out)LeaveOneOut_effectsize <- list()for (i in 1:length(levels(dat$leave_out))) { temp_dat <- dat %>% filter(leave_out != levels(dat$leave_out)[i]) VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5) LeaveOneOut_effectsize[[i]] <- rma.mv(yi = yi, V = VCV_leaveout, random = list( ~ 1 | effectID, ~ 1 | paperID, ~ 1 | species_cleaned), test = "t", sparse = TRUE, data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])}# # writing function for extracting est, ci.lb, and ci.ub from all models# est.func <- function(model) {# df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)# return(df)# }# using dplyr to form data frameMA_oneout1 <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>% bind_rows %>% mutate(left_out = levels(dat$leave_out))# telling ggplot to stop reordering factorsMA_oneout1$left_out <- as.factor(MA_oneout1$left_out)MA_oneout1$left_out <- factor(MA_oneout1$left_out, levels = MA_oneout1$left_out)# saving the runssaveRDS(MA_oneout1, here("Rdata", "MA_oneout1.RDS"))``````{r}#| fig-height: 16MA_oneout1 <-readRDS(here("Rdata", "MA_oneout1.RDS"))# plottingleaveoneout1 <-ggplot(MA_oneout1) +geom_hline(yintercept =0, lty =2, lwd =1) +geom_hline(yintercept = mod1$ci.lb, lty =3, lwd =0.75, colour ="black") +geom_hline(yintercept = mod1$b, lty =1, lwd =0.75, colour ="black") +geom_hline(yintercept = mod1$ci.ub,lty =3, lwd =0.75, colour ="black") +geom_pointrange(aes(x = left_out, y = est,ymin = lower, ymax = upper)) +xlab("Species left out") +ylab("Zr (effect size), 95% CI") +coord_flip() +theme(panel.grid.minor =element_blank()) +theme_bw() +theme(panel.grid.major =element_blank()) +theme(panel.grid.minor.x =element_blank()) +theme(axis.text.y =element_text(size =6))``````{r, fig.width=20, fig.height=30}#| fig-cap: "**Figure S31.** Leave one species out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals."leaveoneout1```<!--  -->## Precision```{r precision}precision_plot<-funnel(mod1, yaxis="seinv", # = "rstudent", xlab = "Standarized residuals", ylab = "Precision (inverse of SE)", ylim = c(0.001, 16), xlim = c(-10,15), digits=c(0,1))``````{r}#| fig-cap: "**Figure S32.** Relationship between standardized residuals and precision (inverse of SE)."precision_plot```:::## R Session Information```{r}# pander for making it look nicersessionInfo() %>%pander()```
